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Abstract 

Background: Assessing heterogeneity in lung images can be an important diagnosis tool. We present a novel and 
objective method for assessing lung damage in a rat model of emphysema. We combined a three-dimensional 
(3D) computer graphics method-octree decomposition-with a geostatistics-based approach for assessing spatial 
relationships-the variogram-to evaluate disease in 3D computed tomography (CT) image volumes. 

Methods: Male, Sprague-Dawley rats were dosed intratracheal^ with saline (control), or with elastase dissolved 
in saline to either the whole lung (for mild, global disease) or a single lobe (for severe, local disease). Gated 3D 
micro-CT images were acquired on the lungs of all rats at end expiration. Images were masked, and octree 
decomposition was performed on the images to reduce the lungs to homogeneous blocks of 2 x 2 x 2, 4 x 4 x 4, 
and 8x8x8 voxels. To focus on lung parenchyma, small blocks were ignored because they primarily defined 
boundaries and vascular features, and the spatial variance between all pairs of the 8x8x8 blocks was calculated 
as the square of the difference of signal intensity. Variograms-graphs of distance vs. variance-were constructed, 
and results of a least-squares-fit were compared. The robustness of the approach was tested on images prepared 
with various filtering protocols. Statistical assessment of the similarity of the three control rats was made with a 
Kruskal-Wallis rank sum test. A Mann-Whitney-Wilcoxon rank sum test was used to measure statistical distinction 
between individuals. For comparison with the variogram results, the coefficient of variation and the emphysema 
index were also calculated for all rats. 

Results: Variogram analysis showed that the control rats were statistically indistinct (p = 0.1 2), but there were 
significant differences between control, mild global disease, and severe local disease groups (p< 0.0001). A 
heterogeneity index was calculated to describe the difference of an individual variogram from the control average. 
This metric also showed clear separation between dose groups. The coefficient of variation and the emphysema 
index, on the other hand, did not separate groups. 

Conclusion: These results suggest the octree decomposition and variogram analysis approach may be a rapid, 
non-subjective, and sensitive imaging-based biomarker for characterizing lung disease. 
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Background 

Emphysema is an obstructive pulmonary disease that re- 
sults in airway expansion and tissue destruction. Early 
and accurate detection of emphysema is important for 
disease management and improved patient outcomes 
[1]. Emphysema typically results in heterogeneous air 
trapping and increased ventilation-perfusion inequality. 
Thus, signatures of tissue heterogeneity or regional air 
trapping may facilitate earlier and/or more accurate 
diagnoses of emphysema [2,3]. 

It has been shown in many studies that CT is a repro- 
ducible and predictive modality for diagnosing and 
assessing emphysema [2,4-7]. The most common way to 
quantify moderate to severe emphysema from CT im- 
ages is to measure the emphysema index: the percentage 
of voxels below a preset Hounsfield Unit (HU) thresh- 
old. However, this index may not detect subtle or early 
onset emphysema [7], or, conversely, it may identify 
"false positive" regions that appear emphysematous in 
asymptomatic healthy young subjects with no smoking 
history [8]. In addition, threshold values vary between 
studies, typically within a range of-900 to-980 HU, with 
factors other than disease influencing CT densitometry 
[2,4,9,10]. 

An approach to assessing heterogeneity that has been 
applied to myriad types of medical images is fractal ana- 
lysis [11]. This approach exploits the scale-independent 
nature of fractals in systems, such as lung vasculature, 
whose variation in form or regularity is thought to be 
similar through different degrees of magnification [12]. 
Typically, a region (or regions) of interest (ROI) is se- 
lected, which is then evaluated for a mean value and/ 
or repeatedly subdivided to relate differences in signal 
intensity across space and across ROI size. If the rela- 
tionship is characterized by a power law, then the ex- 
ponent, or fractal dimension, is taken as an indication 
of complexity, texture, or heterogeneity [13-15]. How- 
ever, this method is based on an a priori assumption 
that power law relationship exists and that a fractal 
dimension can be determined. Furthermore, because 
this approach requires ROI selection, it is often sub- 
jective and can result in omission of large sections of 
the lung [13,16]. 

A promising new approach, detailed in recent work by 
Subramaniam et al. [17], demonstrates the use of ana- 
lysis of CT image slices using a quadtree decomposition. 
This approach iteratively divides a 2D image into in- 
creasingly small squares. Division occurs whenever the 
range of intensities within the quadrant are above a 
pre-defined threshold. In this way, regions with homo- 
genous intensity are characterized by larger squares. Their 
work measures heterogeneity by counting the number 
of squares-per-area in the lung or in a local region of 
the lung. 



We propose to extend the quadtree decomposition to 
3D by using octrees. In this way, the entire lung CT 
image volume is used. An octree iteratively divides a 
cubic volume into eight evenly sized cubes, or octants. 
The octree method is well-developed and is used, for ex- 
ample, in 3D mesh generation and computer graphics 
[18,19], medical image registration [20], and finite elem- 
ent meshing in CT images [21], The advantages of the 
octree here is that it facilitates a speed-up in 3D lung tis- 
sue analysis by non- subjectively subdividing the lungs 
into homogeneous regions, thereby allowing for focus 
on parenchyma by reducing partial volume averaging 
and eliminating edges [22]. 

To incorporate spatial information across the lung, we 
propose to couple the octree image decomposition ap- 
proach with variogram analysis. Variograms are a spatio- 
statistical approach for measuring spatial variability by 
comparing sample value variances to the distance of 
separation [23,24]. Variograms are well-established in 
geostatistics, but have also recently attracted interest in 
biomedical applications, such as characterizing magnetic 
resonance images of white matter [25]. In this pilot study 
using our approach, we are able to significantly differenti- 
ate from each other the three groups of elastase-dosed 
rats: control, distributed mild emphysematous disease, 
and region-specific severe disease. 

Methods 

Disease model 

An elastase-induced model of emphysema was used in 
this study. Nine male Sprague-Dawley rats with an aver- 
age weight of212±ll g were orally intubated and dosed 
intratracheally with: 250 U/kg elastase dissolved in 
200 [iL saline to the whole lung (n = 3), or 50 U/kg elas- 
tase in 200 [iL saline to a single lobe (n = 3), or 200 \iL 
saline as a control (n = 3). Dosing levels were based on 
our previous work in which emphysematous changes 
were detected using 3 He diffusion MRI and histology 
[26]. All animal use was approved by the Institutional 
Animal Care and Use Committee at Pacific Northwest 
National Laboratory. 

CT imaging 

Three weeks following dosing, the rats were imaged 
using micro-CT. At this time, rats weighed 357 ± 10 g. 
The imaging procedure is described in detail in [27]. 
Briefly, rats were anesthetized, intubated, and mechanic- 
ally ventilated at 1 Hz with 40% inhale and 60% exhale 
durations. Peak inhalation pressure was ~8 cmH 2 0, and 
no peak end expiratory pressure was used so that images 
could be acquired at functional residual capacity. Anes- 
thesia was maintained by providing 3-4% isoflurane in 
air (30% 0 2 , balance N 2 ). Sigh breaths were delivered 
periodically to maintain lung recruitment. A respiratory- 



Jacob and Carson BMC Medical Imaging 2014, 14:1 
http://www.biomedcentral.eom/1471-2342/14/1 



Page 3 of 1 1 



gated GE eXplore 120 micro- CT scanner was employed 
with the following settings: 100 kV peak voltage, 50 mA 
tube current, 16 ms exposure time, and 360 projections 
with 1 degree angular steps. Images were reconstructed 
with supplied software to 150 \im isotropic resolution. 
Total imaging time was about 90 minutes due to the 
collection of multiple images throughout the breathing 
cycle; however, only the images acquired at full exhal- 
ation were used for the analyses herein. Post-mortem 
histology and tissue analysis were not possible due to in 
situ lung casting for other research purposes (i.e. supply- 
ing airway tree geometries for computational fluid dy- 
namics models [28,29]). 

Image preparation 

A lung mask image was semi-automatically generated 
from the above-mentioned reconstructed images using 
Image J [30] and the Image J 3D Toolkit plug-in [31]. 
Starting from a seed-point inside the lung, a 3D con- 
nected threshold was applied with a threshold value em- 
pirically selected to exclude major vasculature and all 
external tissue, so that only lung tissue would be in- 
cluded in analysis. Mask boundaries were smoothed 
using the Region Dilate and Region Erode functions in 
succession. Applying the mask to the reconstructed im- 
ages, all background was thus assigned an intensity value 
of 0, and all unmasked lung tissue retained original HU 
values. A five pixel diameter 3D median filter was then 
applied in order to reduce noise while contributing 
minimal blurring. Finally, the image canvas size was 
increased to 512 x 512 x 512 by zero-filling in each 
direction. We note that zero-filling to a power of 2 
served to conveniently restrict all octants in the octree de- 
composition to isotropic cubes. 

Octree decomposition 

An automated octree decomposition was performed by 
iteratively subdividing an image, with each division pro- 
ducing eight evenly sized octants (e.g., the initial 512 x 
512 x 512 image was decomposed into eight 256 x 256 x 
256 octant regions, and so on; see Figure 1) [32]. After 
each division, the maximum, minimum, mean, and stan- 
dard deviation of the signal intensity were calculated for 
each octant region. Then, iterative octant subdivision of 
a region occurred if either the standard deviation of 




that region exceeded a user-defined threshold, or if the 
region contained at least one voxel-but not all voxels- 
with intensity equal to the background signal (i.e. 0). 
However, regions reaching a minimum size of 2x2x2 
were not subdivided. At the end of the decomposition, 
all boxes with a mean HU value * 0 and dimensions 
greater than 2x2x2 represented portions of lung with 
relatively uniform HU values. Following decomposition, 
results were binned according to HU and box size for 
histogram display. Boxes containing at least one 0 were 
not included in the binning. Octree decomposition was 
implemented in Python and executed on a MacPro 
model 3.1. We note that we have modified our Python 
code to integrate ImageJ to minimize user interaction, 
although that was done subsequent to this work. The 
only interactive user-initiation required is the manual se- 
lection of a seed-point for the 3D connected threshold. 

There are multiple approaches to selecting the afore- 
mentioned user-defined threshold criteria. For example, 
in [17] boxes were subdivided if the HU range in the 
box exceeded 10% of the range of values in the entire 
image. Because our rats had different ranges of HU 
values depending on disease severity, we chose to use 
the standard deviation as the threshold criterion, as it 
does not require normalization. We describe how the 
threshold range was determined in Section Threshold 
range. 

Variogram analysis 

The purpose of using variograms in this context was to 
determine the spatial relationships between regions of 
the lung that are relatively homogeneous. The octree de- 
composition isolated regions of relatively high homo- 
geneity; therefore, each resulting box was treated as a 
homogeneous "voxel" using the mean value of the box 
as the voxel intensity. Furthermore, it was important to 
use octree boxes that were all the same size. The largest 
box sizes that resulted from the octree decompositions 
were mostly 8x8x8 with a few 16x 16x 16 boxes. In 
order to ensure a sufficient number of boxes that were 
distributed throughout the lung and to enhance sta- 
tistical power [25], each 16x16x16 box was further 
broken down into 8x8x8 boxes, and all 8x8x8 
boxes were used in the variogram analysis. This ap- 
proach assumed that emphysematous disease was ma- 
nifest on a scale greater than an 8x8x8 box, which 
for these images is 1.2 mm on a side-this is consis- 
tent with what we observed in this disease model in 
previous work [26]. Importantly, excluding the smaller 
boxes largely eliminated vasculature, conducting air- 
ways, boundaries, and other features that had high 
spatial variability [22]. 

The calculation of the variance (or semi-variance) y(d) 
of the differences in signal intensity / is described in 
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detail in many geostatistics texts, such as [23,24], and is 
given by: 

y(d) = -^Er^)-^-^)] 2 (i) 

where d is the distance between voxel centroids, x is the 
voxel location in 3D, and N(d) is the number of voxel 
pairs for a given d. y(d) is calculated for all values of d. 
For an image with n voxels, the total number of voxel 
pairs (for all existing values of d) is given by 

N = n(n-l)/2. (2) 

We then calculated the average semi-variance at each 
distance. Results were plotted on a distance vs. variance 
graph, or variogram, which graphically represents the 
spatial dissimilarity within the image. 

The purpose of the variograms was to characterize 
spatial relationships between boxes of relatively high 
homogeneity that were dispersed throughout the lung. 
However, as described by [25] in work with brain im- 
ages, increased distances between voxels, or octants, lead 
to decreased likelihood that the voxels are related in any 
way. This is especially true in the lung, which is com- 
posed of largely independent lobes (five in the rat). 
Moreover, the spatial relationship between parenchymal 
signal intensity in different lobes is, in many contexts, 
most accurately described by a unique path up and back 
down the airway tree, potentially covering 20 or more 
orders of branching [33]. The spatial distances are analo- 
gous to the "regionalized variable" in the geological con- 
text of mineral distribution. Indeed, the regionalized 
variable, or the maximum distance d max over which the 
variogram is expected to be reliable, is estimated to be 
one half of the diameter of a region [25], while beyond 
d max any relationships between distance and variance 
are expected to be random. Therefore, we limited our 
variogram analysis to a region that represented half of 
the average dimension of the left lobe, which is typically 
the largest lobe in the rat. This assumes that the disease 
varies slowly compared to the box size but rapidly com- 
pared to the lobe size. Using the 3D images and a lung 
cast, d max was measured to be ~ 8 mm (53 pixels) in the 
healthy rats, which was rounded up to the equivalent of 
seven face-bordering 8x8x8 boxes, or 56 pixels. 

Octree decomposition tests 

We performed the following evaluations of the appro- 
priateness and robustness of the octree and variogram 
approach. 

Threshold range 

The octree decomposition compared the standard devi- 
ation of the lung signal within a box to a threshold range 



to determine whether the box should be subdivided. The 
optimal threshold range was determined semi-empirically. 
As a starting point, we calculated the typical standard de- 
viation of the lung tissue of the three control rats. This 
was accomplished by producing a histogram of each 
masked image and then fitting the main lung peak to a 
Gaussian curve. The mean of the control rats' standard 
deviations (a) was set as the initial threshold range. We 
then tested the effects of thresholds that were a/3, 2o73, 
4o73, and 2a by rerunning the decomposition on the same 
images. We evaluated the effectiveness of each threshold 
range at grouping the control rats and at distinguishing 
the full-lung-dose rats from the control group in the 
variograms. 

Image filtering 

We examined the effects of image filtering by generating 
variograms from unfiltered images, images filtered with 
a five pixel diameter 3D median filter, images filtered 
with a 3D Gaussian blur of radius = 2, and images fil- 
tered with a 3D Gaussian blur of radius = 4. Filters were 
applied in ImageJ prior to applying the mask. The num- 
ber of octree boxes and differences in the variograms 
were compared. Adjustments to the threshold range 
were made for each filter test based on the standard de- 
viation of the filtered image. 

Image translation 

Rat positioning during imaging can vary from animal to 
animal, and the octree decomposition should be inde- 
pendent of this. Therefore, we tested the effects of image 
translation on the variogram results. For the variogram 
analysis, the largest unit box retained after decompo- 
sition was 8x8x8; therefore, a shift in the image by 8 
pixels in any (or all) dimensions should result in no 
change to the resulting variogram, since the original 
image had isotropic dimensions of 2 n . Conversely, ma- 
ximum change would occur from a four-pixel shift. To 
test the significance of the effects of translation, the 
image of one rat was shifted in x, y, and z by four pixels 
and by eight pixels using ImageJ, and variogram results 
were compared to those of the original image. 

Image rotation 

Animal positioning can also have an apparent effect on 
image rotation, and octree decomposition should be 
insensitive to this. Using the same rat image as used 
in the translation test, we imitated an arbitrary 3D 
rotation by rotating the image n/4 radians about the 
x, y, and z axes using ImageJ. Results were compared 
to the original image and to the results of the transla- 
tion test. 
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Image downsampling 

To demonstrate the increased value of the combined 
octree and threshold approach for establishing the 8 x 
8x8 boxes versus simply downsampling the image, we 
generated variograms using standard image downsam- 
pling as a comparison. For this, we applied a bilinear 
downsampling to the 512x512x512 image, creating a 
64 x 64 x 64 image-each voxel the equivalent size of an 
8x8x8 octree box. Then we generated variograms utiliz- 
ing the intensity value of every voxel in the 64 x 64 x 64 
image (excluding the background voxels, again defined as 
those voxels with intensities equal to zero). 

Emphysema index, coefficient of variation, and 
heterogeneity score 

The percentage of lung below a HU threshold value, or 
emphysema index, was calculated for each rat to com- 
pare this conventional measurement of disease severity 
with the variogram results. Because there is no estab- 
lished HU threshold level in rat models of emphysema, 
we chose to count the percentage of voxels with HU 
values below two standard deviations from the control- 
group mean. This level was determined to be-717 HU. 
Calculations were made on the same masked images used 
for octree decomposition. 

The CoV was calculated for each rat by fitting a histo- 
gram of the masked lung images to a Gaussian curve 
and taking the ratio of the standard deviation to the 
mean. 

For comparison, we defined a new metric, the hetero- 
geneity score (A), as the average difference between the 
mean variance of the control group and the spatial var- 
iance of each individual rat, in the range d < d max . 

Statistical analysis 

Comparisons of variograms of the three control group 
rats were made using the Kruskal-Wallis (KW) rank sum 
test with a null hypothesis of a = 0.05. This indicated 



whether there were significant differences within the 
group. For pairwise comparisons, a Mann-Whitney- 
Wilcoxon (MWW) rank sum test was employed, also 
with a null hypothesis of a = 0.05. 

Results 

CT images 

Figure 2 shows a representative unfiltered coronal slice 
from a rat in each dose group: panel A is a control rat, 
panel B is a full-lung-dose rat, and panel C is a partial- 
lung- dose rat (the dose was delivered to the distal por- 
tion of the left lobe). The mean HU value (± standard 
deviation) of the control group was-544 ± 58 HU and 
the mean of the full-lung dose group was-594 ± 38 HU. 
In spite of marginally lower HU values that would be 
expected from emphysematous disease, an analysis of 
variance showed that the two groups were indeed not 
statistically different (p = 0.26). On the other hand, the 
partial-lung-dose rats showed a bimodal distribution of 
HU values because the diseased regions of the lung were 
distinct from the healthy regions. An example of this 
can be seen in the lower portion of the left lobe of the 
partial-lung-dosed rat (panel C) where the signal in- 
tensity is substantially lower than the rest of the lung 
(by -200 HU) -indicative of severe tissue destruction 
and air trapping, which are characteristics of emphyse- 
matous lungs [34]. Based on the overall HU measure- 
ments, we presume that the full-lung- dose rats developed 
a mild-and difficult to distinguish-emphysematous dis- 
ease while the single-lobe-dose rats developed a more se- 
vere, albeit localized, disease. 

Octree decomposition tests 
Threshold 

The starting octree decomposition threshold level deter- 
mined from the standard deviations of the control rats 
(Section Threshold range) was 60 ± 1 HU; therefore, 60 
HU was used as the initial threshold level. Thus, the 



A A A 

Figure 2 Unfiltered coronal slices from 3D images of A) a control rat, B) a full-lung-dose rat, and C) a single-lobe-dose rat. Only subtle 
differences between the full-lung-dose and control rats are evident; however, the distal region of the left lobe of the single-lobe-dose rat has 
considerably lower signal intensity, indicative of tissue destruction and/or air trapping characteristic of emphysematous disease. 
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other threshold levels tested were 20, 40, 80, and 120 
HU. Octree decompositions were performed using the 
different threshold levels on the three control rats, and 
variograms were compared using a KW rank sum test. 
The threshold level that showed the least difference 
among the control rats was 40 HU (p = 0.12), with 60 
HU also showing no significant differences (p = 0.09). 
The 20, 80, and 120 HU thresholds did have significant 
differences among the controls (p = 0.02, p < 0.0001, and 
p < 0.0001, respectively). Based on this result, we chose 
to use the threshold level of 40 HU for all subsequent 
octree decompositions (unless otherwise noted). 



Image filtering 

Results of applying the different filters showed that the 
variograms from the unfiltered image and the image 
with the median filter were indistinct in an MWW test 
(p = 0.62). However, the relatively noisy unfiltered image 
resulted in about 20% fewer 8x8x8 boxes than the fil- 
tered image. This was in spite of a higher threshold that 
was used, 48 HU, based on the unfiltered images stan- 
dard deviation. The radius = 2 and radius = 4 Gaussian 
filter results also did not differ significantly from the 
unfiltered image (p = 0.41 and p = 0.48, respectively) or 
from the median-filtered image (p = 0.71 and p = 0.84, 
respectively). A KW test showed that the radius = 2 
filter resulted in control group variograms that were 
not distinct (p = 0.16), but the radius = 4 variograms 
were (p = 0.0001). However, because of the blurring 
caused by Gaussian filters, the number of 8x8x8 
boxes that resulted from the Gaussian-filtered images 
was about 30% higher than the median filtered image, 
presumably because vascular structures and edges 
were not well preserved. Threshold levels used for 
the radius = 2 and radius = 4 images were 37 HU and 
39 HU, respectively, based on the post-filtering stan- 
dard deviations. 



Effect of downsampling 

Figure 3 shows an example of an original image (panel 
A), the image with the five voxel diameter 3D median 
filter (panel B), and the same image downsampled to 1/8 
resolution, or 64 x 64 x 64 pixels (panel C). The vario- 
grams made from the downsampled images of the three 
control rats were statistically different in a KW rank 
sum test (p < 0.0001). On the other hand, the vario- 
grams made from only the 8x8x8 boxes that came 
out of the octree decomposition were not (p = 0.12; 
see Section Threshold). 

Translation and rotation 

Results of image translation and rotation were compared 
for one rat. For the eight-pixel translation, variogram re- 
sults were exactly identical to that of the original image, 
as expected. The four-pixel shift was not statistically dis- 
tinct from the original image as confirmed by a MWW 
test (p = 0.96). In addition, the result of the rotation 
showed no significant change from the original (p = 0.56) 
or from the four-pixel shift (p = 0.59). Thus, we confirm 
that shifts or rotations to the image (i.e. alternative posi- 
tions of the lung during imaging) do not result in signifi- 
cant changes to the resulting variograms. 

Octree decomposition 

Figure 4 shows the results of the octree decomposition 
on a control rat and on one with severe disease in the 
lower left lobe (see Figure 2C). Column A shows the 2 x 
2x2 boxes, and column B shows the 4x4x4 boxes. 
These box sizes largely define the fine structures and 
edge details, including the conducting airways and vas- 
culature. For this reason, we ignored the 2x2x2 and 
4x4x4 boxes for the variogram analysis. Column C 
shows the 8x8x8 boxes and their relatively uniform 
distribution throughout the lung. Conversely, the less 
uniform distribution of the considerably fewer 16 x 
16 x 16 boxes-the largest that resulted from octree 
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Figure 3 A masked axial slice from an original CT image with examples of different filtering. A: Unfiltered CT image. B: 3D median filter 
(5 pixel diameter). C: Downsampled by 8x. 
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decomposition-is shown in Column D. We point out, for 
example, that the lower left lobe of the treated rat (column 
D, bottom) has an approximately 2.5 x higher density of 
16x16x16 boxes than that of the control rat, indicating 
prevalence of localized homogeneous CT signal intensity. 

Histograms of the octree decomposition for three rats 
are shown in Figure 5. The marker size corresponds to 
the box size. In the control and mild disease rats (panels 
A and B), the different sized boxes, other than the 2 x 
2 x 2s, are approximately Gaussian distributed about the 
mean HU value of the lung, whereas the single-lobe-dose 
rat (panel C) shows a bimodal distribution, indicative of at 
least one large region of the lung with considerably lower 
HU values. These histograms show what percentage of the 
lung at each HU intensity is defined by the different box 
sizes, but they do not convey any information about the 
spatial relationships of the boxes. 

Variograms 

In order to visualize the spatial relationships between 
the 8x8x8 boxes (and decomposed 16 x 16 x 16 boxes), 
variograms were constructed. Figure 6 shows the average 
variograms from the three different dose groups. The 
dashed line in Figure 6 denotes the range of d max (see 
Section Variogram analysis). For d < d max a KW test 
showed that the control rats were statistically similar 
(p = 0.12), and a MWW test showed that the dose 



groups were each statistically distinct from the control 
group (p < 0.0001) and from one another (p < 0.0001). 

Emphysema index, coefficient of variation, and 
heterogeneity score 

Results of the emphysema index are shown in the 
Figure 7A. Although the number of rats is too small 
to reliably calculate sensitivity and specificity, the graph 
shows that there is considerable overlap between subjects 
in each dose group, likely indicating poor sensitivity and 
specificity. The CoV for each rat is shown in Figure 7B. 
Similarly, there is no distinction between the control and 
full-lung dose groups. The part-lung dose group had con- 
siderably higher CoV values, because the standard devi- 
ation was enlarged due to a bimodal (and non-Gaussian) 
distribution of HU values. 

For comparison, we calculated the heterogeneity 
score A of each rat s variance from that of the control 
group average; see Figure 7C. There is no overlap of A be- 
tween dose groups (i.e. the highest A of the control group 
is lower than the lowest A of the full-lung group). This is 
not observed in Figure 7A or B, suggesting better sensitiv- 
ity and specificity for A than the other metrics. 

Discussion 

This work combined two different data analysis tools, 
octree decomposition and variograms, to study tissue 
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Figure 5 Octree decomposition histograms showing the 
relationship between box size, HU, and percentage of the lung 
by each box size at each HU level for three representative rats. 

A: Control (mean = -599 HU), B: Full-lung dose (mean = -634 HU), 
C: Single-lobe dose (mean of dosed region = -801 HU, mean of 
undosed region = -530 HU). 



heterogeneity in lung disease. We showed that this 
merged approach was better able to differentiate rats 
with mild emphysematous disease from the healthy 
control group than methods that relied on absolute 
HU values. The main criterion for octree decomposition 
was based on the standard deviation of HU values within 
an octree box. An advantage to this approach is that it 
avoids thresholding according to HU values, although so- 
phisticated thresholding algorithms may be useful [35-37]; 
rather, it focuses only on heterogeneity-based signatures 
that may characterize disease [2]. We propose that a het- 
erogeneity score A, the average distance of a rats variance 
from that of the control group average, may be useful to 
classify disease severity. Furthermore, to visualize the re- 
gions of the lung with the greatest heterogeneity, one 
could determine which boxes had the highest semi- 
variance within d max and map them back to the ori- 
ginal image. This would provide 3D information about the 
spatial distribution of lung tissue heterogeneity and, po- 
tentially, disease distribution. 

Another approach to a disease metric might be that of 
fitting the data to an established variogram model, most 
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Figure 6 Variograms averaged for each dose group. The dashed 
line shows the extent of d maxi the characteristic lung distance to 
which variograms were analyzed for heterogeneous variations. The 
lines get "noisier" at higher distances, indicative of increased spatial 
variation even within the each group. 
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Figure 7 Graphs comparing the different metrics for detecting 
emphysema. A: Percentage of image pixels with HU values below 
two standard deviations from the mean of the control group. 
B: Coefficient of variation (CoV). C: The average distance of each 
individual rat's variogram from the mean of the control group. The 
x-axis labels are animal numbers: 100s are controls, 200s are full- 
lung-dose, and 300s are single-lobe-dose. 301 and 302 received the 
dose in the right caudal lobe, and 303 received the dose in the 
distal region of the left lobe (Figure 2C). 



of which describe an asymptotic rise in variance (i.e. 
variance becomes independent of distance indicating 
that spatial relationships become random) [23]. This is 
seen to some degree in our data (see Figure 6). However, 
within the range of d < d max , we found that a power law 
model generally fit our data better than other models. 
This model typically describes fractal behavior [38]. Our 
initial investigations into this model showed potential 
promise at using fit parameters to distinguish dose 
groups; however, results lacked statistical significance, 
and the full-lung dose group tended to not fit this model 
as well as the other two groups. 

There was no single variogram model that satisfactor- 
ily fit the all the data over the entire range of distances, 
because the complex geometries found in the lung result 
in some problems for variograms. In particular, the dir- 
ect linear path between two regions of the lung sepa- 
rated by large distances often crosses non-lung tissue, 
such as the heart, which are essentially treated as holes 
or voids in the geometry. Furthermore, neighboring lobes 
generally do not interact physiologically except through 
the vascular and airway trees, which may only connect re- 
gions through many orders of branching. As pointed out 
by Keil et al. [25], one could go to extraordinary measures 
to take into account structural distances (the physio- 
logical distance at which different regions interact) versus 
Euclidian (straight line) distances used herein. To limit the 
problem in this study, we constrained the distance of var- 
iogram analysis to approximately half the characteristic 
diameter of the largest lobe. However, to better under- 
stand the inter-and intra-lobe variance relationships, the 
lung could be segmented into lobes (if the image is of 
sufficient resolution to discern lobar boundaries), and the 
decomposition/variogram process repeated on the seg- 
mented images. Though, by ignoring inter-lobar variances, 
this approach would likely not capture information about 
disease that was confined to a single lobe, particularly if 
the entire lobe was affected homogeneously. 

We employed octree decomposition prior to genera- 
ting the variograms. In doing this, we assumed that the 
emphysematous disease is generally slowly varying over 
space, and that the disease causes changes to homoge- 
neity on the order of or less than the lobar length scale 
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but greater than the octree box length scale [4,34]. This 
is consistent with what we previously observed using 
3 He MRI in the same disease model [26]. Without this 
assumption, variograms would have to be made directly 
from the raw images. This is possible but impractical, 
because the semi-variance computation time (and re- 
sulting file size) is proportional to the number of voxel 
pairs, which rises approximately as n 2 (see Eq. 2). We 
measured the variogram computation time for 1078 
octree-decomposed boxes from one rat to be 6.6 se- 
conds (on a MacPro model 3.1), and we verified ex- 
perimentally that the computation time indeed rose 
in proportion to the square of the number of voxels. 
Therefore, to create a variogram on the entire masked 
3D image, which consisted of 1.48 x 10 6 voxels (of 
lung tissue only), it would take us ~ 1.25 x 10 7 seconds, or 
about 5 months-with a resulting file size on the order of 
20 GB. Therefore, octree decomposition dramatically 
reduces the computation time while focusing non- 
subjectively on regions of the lung that are of greatest 
interest. We note that the octree decomposition itself 
was performed in ~2 minutes. 

Image noise can confound octree decomposition and 
affect resulting variograms. Results of image filtering tests 
indicated that noise reduction using an edge-preserving 
filter resulted in more 8x8x8 octree blocks without sig- 
nificantly affecting variogram results. An alternative to 
octree decomposition is downsampling, which is a quick 
and straightforward approach to reducing image noise 
and size. However, we verified that the octree decompos- 
ition approach performed much better at separating dose 
groups than simply downsampling the image and then cal- 
culating the semi-variance using every voxel. The octree 
decomposition assures that only the homogeneous regions 
of the lungs are singled out for comparison, whereas 
downsampling the image blurs together proximal voxels, 
including vasculature, airways, and lung boundaries 
irrespective of signal intensity or tissue type. Thus, the 
downsampling approach apparently causes a loss of infor- 
mation. The radius = 4 Gaussian filter had a result similar 
to downsampling. 

One limitation of this pilot study was the small num- 
ber of animals in each group, which did not allow the 
statistical evaluation of specificity and sensitivity. There- 
fore, follow-on work will be required to validate these 
results and establish specificity and sensitivity [2]. This 
might be accomplished in conjunction with pulmonary 
function tests, conventional morphometric measure- 
ments [39], and histological techniques particularly sen- 
sitive to early emphysematous changes [40]. Additional 
future work should evaluate the performance of this 
method on clinical CT images as well as test the effec- 
tiveness for distinguishing different diseases and disease 
models. 



Conclusion 

Results of this pre-clinical study of elastase-treated rats 
suggests that automated octree decomposition and vario- 
gram analysis based on image heterogeneity may provide 
a non-objective and sensitive metric for characterizing 
emphysematous lung disease, even in early disease stages. 
The method outperformed conventional approaches that 
utilize thresholding and absolute HU values. This ap- 
proach may be applicable to human datasets and other 
diseases. 
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